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Theoretical continuum models that describe the formation of patterns on surfaces of targets under- 
going ion-beam sputtering, are based on Sigmund's formula, which describes the spatial distribution 
of the energy deposited by the ion. For small angles of incidence and amorphous or polycrystalline 
materials, this description seems to be suitable, and leads to the classic BH morphological theory 
[R. M. Bradley and J. M. E. Harper, J. Vac. Sci. Technol. A 6, 2390 (1988)]. Here we study the 
sputtering of Cu crystals by means of numerical simulations under the binary-collision approxi- 
mation. We observe significant deviations from Sigmund's energy distribution. In particular, the 
distribution that best fits our simulations has a minimum near the position where the ion penetrates 
the surface, and the decay of energy deposition with distance to ion trajectory is exponential rather 
than Gaussian. We provide a modified continuum theory which takes these effects into account and 
explores the implications of the modified energy distribution for the surface morphology. In marked 
contrast with BH's theory, the dependence of the sputtering yield with the angle of incidence is 
non-monotonous, with a maximum for non-grazing incidence angles. 

PACS numbers: 05.10.-a, 68.35.-p, 79.20.-m 



INTRODUCTION 



Ion bombardment of solids often gives rise to char- 
acteristic surface topographies, which evolve under sta- 
tionary and homogeneous ion fluxes. Besides kinetic 
roughening, wavelike ripple structures may occur. Such 
height modulations on the submicron scale have been 
observed for crystalline semiconductors^ as well as for 
crystalline metals 3 ^ and some amorphous 5 and polycrys- 
talline materials, see a recent review in Ref. |(J Ac- 
cording to continuum theories, which are based on the 
work of Bradley and Harper (BH), 7 the periodic pat- 
terns emerge from a competition between a roughen- 
ing curvature instability due to characteristics of the 
spreading of ion energy, and simultaneous smoothing pro- 
cesses due to surface diffusion*^ Although this mech- 
anism seems to be quite universal, there are material- 
specific differences in the evolution of surface topogra- 
phies. For non-metallic substrates, for example, one 
usually needs off-normal incidence of ion flux to pro- 
duce ripples, which change their orientation with the 
incidence an gl e j2i!i£ii2iiiii2iiiii4 whereas ripples are ob- 
served on metallic substrates even at normal incidence, 
and the orientation of ripples changes with substrate 
temperatureAi 5 ^^ Furthermore, the smoothing mecha- 
nism of surface diffusion is not well understood yet. In 
previous simulations^ we have found that the emerging 
patterns depend crucially on the diffusion mechanisms 
applied. In particular the long-time behavior, which is 
governed in the continuum theory by non-linear terms, 
depends even qualitatively on the surface diffusion mech- 
anism. Given that the surface topographies resulting 



from different mechanisms of surface diffusion have been 
studied by simulations elsewhere^ in the present work 
we will focus on specificities due to the energy deposition 
process. 

Continuum theories for the surface morphology of the 
target usually assume that the kinetic energy of an ion 
hitting a solid surface spreads in the bulk and produces 
a Gaussian density of deposited energy 



■'■ - - <i 



e s (r) = N s ee 



e 

21-1 



(1) 



where r = (x, y, z) is a point within the target, ions are 
falling along the z axis and penetrate an average distance 
a within the solid, e is the average kinetic energy carried 
by each ion, and the values of a, (3 describe the spread- 
ing of the energy, they are of the same order of magni- 
tude as a. The Gaussian form is based on the work 
of SigmundjA^ who considered a polycrystalline or amor- 
phous target and analyzed the kinetic transport theory 
of the sputtering process. He found that in the elastic 
collision regime at energies where electronic stopping is 
not dominating, the deposited energy can be approxi- 
mated by a Gaussian near its maximum. The quality of 
the approximation is reasonable, if mass differences be- 
tween substrate and ion are not too large. Obviously, 
the Gaussian form is not universal and consequences of 
deviations from the Gaussian form within the BH model 
have not been studied yet. In particular, although the 
observations of ripples on single crystalline metals 315 - 16 
are qualitatively described by the BH model, the latter is 
strictly a theory for amorphous materials, and thus there 
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is a need to justify theoretically the emergence of such 
type of patterns onto this other class of substrates. 

Obtaining more detailed information about the de- 
posited energy from simulations has become straightfor- 
ward by now, as there are many well calibrated, efficient 
simulation methods for ion impact availablei^fi^i'SSi2iS4 
In the present work, we use simulations based on the bi- 
nary collision approximatiori22i22i2^ and consider a metal- 
lic material (Cu), for which we generate statistical ensem- 
bles of collision cascades emerging from single ion impact 
events on plane surfaces. We analyze the data in terms 
of deposited energies as well as ejected particles. We do 
not claim to perform the best state of the art simula- 
tion of ion impact on Cu (for example, we only consider 
a very rough model of surface binding forces). Rather, 
we aim at more generic results, which are of relevance to 
the theory of surface evolution. Our simulations provide 
an average density of deposited energy, which is quite 
different from the Gaussian form in Eq. . We further- 
more consider the fluctuations around this average and 
find strong, intrinsic noise. In the subsequent part of our 
work we investigate the consequences of the simulation 
results for the continuum theory of pattern formation by 
ion-beam sputtering. We obtain that the modified en- 
ergy distribution obtained in the numerical simulations 
induces a sputtering yield that overcomes some of the 
shortcomings (when comparing with experiments) of the 
analogous result within BH's theory. Moreover, we re- 
cover the production of the ripple instability, and the de- 
pendence of the pattern features with phenomenological 
parameters similar to BH theory, thus providing a the- 
oretical framework within which observations of ripples 
on metals can be naturally accomodatedAi^i^ 

II. OBSERVABLES OF CASCADE STATISTICS 

In this section we want to relate observables of our 
simulations to the phase space density 

g(v,p,z,t\p o ,z = 0,v ,t = 0), (2) 

where p = (x, y). This function is the basic quantity un- 
derlying the kinetic theory of collision cascades and also 
introduces the quantities which are used in the construc- 
tion of a continuum theory of surface pattern formation 
by ion bombardment. Function g is the average density 
of cascade particles in 6-dimensional (v,p, z) space at 
time t, under the assumption that one ion has hit the 
surface at po and at t = with velocity Vq- As we will 
only treat identical initial conditions with po — and 
■Wo = — |^/2eo/m|e 2 , we will use the abbreviated nota- 
tion g{p, z,v,t) and drop the explicit dependence on po 
and Vo- The average has to be taken over an ensemble of 
targets, which differ by random, thermal displacements 
of atoms. 

To define our simulation observables in terms of 
the phase space density, first note that g(p, z,v,t)(v ■ 
da) cPvdt is the number of particles, which penetrate a 



surface element da situated at position (p,z), with ve- 
locity v during the time interval dt. 

The phase space density and the corresponding current 
density may as well be considered as functions of position, 
energy and direction of velocity using v = yj2t/mv, v 
being the unit vector in the direction of v, so that vd 3 v = 
(2/m 2 )ev dedv. The current density j(e, p, z, t)de of cas- 
cade particles of energy near e is given by 

j(e, p, z,t)de = — -de [ dv evg(e,v, p, z,t). (3) 
m z J 

If we restrict the integration over v to directions with 
v-da > 0, we only count particles, which cross the surface 
element in the direction of its normal. This variant will 
be denoted by j + . 

From Eq. © it is obvious that the time integral 

/>OC 

h 2d (e, p, z) dxdy = / dt e z ■ j(e, p, z,t) dxdy (4) 
Jo 

equals the total average number of particles per energy 
at energy e of a single collision cascade, which penetrate 
the surface element dxdye z located at (p, z). 

Note that h 2 d is a surface density and the quantities 

/>oo 

n 2d (p,z) = / de h 2d {e,p,z) (5) 
Jo 

and 

I* oo 

e-2d{p, z) = / de eh 2d {e,p,z) (6) 
Jo 

give the average number of particles per area and average 
energy per area transported to the (ccy)-plane at z by the 
collision cascade. For all these quantities the correspond- 
ing variants h^ d , n^ d , e^ d only take into account particles 
moving in outward direction +e z . 

The particles arriving at the z = plane (which con- 
stitutes the surface of the material) with velocities in 
outward direction will leave the bulk if they overcome 
the surface binding forces. We will use a simple spherical 
barrier model of surface binding with barrier height U. 
This implies that all particles arriving at the surface with 
kinetic energy e > U will be sputtered off. The surface 
density njj of these particles is therefore given by Eq. 

with the lower boundary of the e integration replaced 
by U. The total sputtering yield is the surface integral 
of this density, Yjj — j d 2 p nu{p). At internal surfaces 
there is no surface binding and thus 

p(e,p,z) = h 2d {e,p,z)/Y u=Q (7) 

becomes the probability density to find a particle with 
energy e crossing the internal surface at location (p, z) 
and Pu(p) — n u{p) /Yjj is the probability density to find 
a particle leaving the bulk at p. These are the quantities 
we will study in the subsequently described simulations. 
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III. BCA SIMULATIONS 

Atomic displacements and particle ejection from a solid 
due to the impact of a single ion with kinetic energy in 
the keV range can be simulated by using the binary col- 
lision approximation 24 (BCA). The basic idea is to sub- 
stitute the detailed particle trajectories by trajectories 
where the particles travel with constant velocity until 
they "hit" onto another particle. Each collision event 
is integrated analytically or numerically, leading to new 
positions and velocities of the particles participating in 
the collision. Hence the full dynamical process is reduced 
to a cascade of collisions. A sample cascade, originating 
from an impact of a 5 keV Cu ion on a Cu lattice with an 
angle of incidence of 60° , is shown in Fig. Although 
this method has its limitations^ it has become a stan- 
dard technique and is used to describe ion implantation 
and sputtering. 



microscopically differing configurations of the solid. 

We considered two orientations of the crystal, (1, 0, 0) 
and (58, 72, 39). The latter was used to suppress effects 
of crystal anisotropy and we found very good agreement 
between the angular averages of ri2d, &id obtained from 
(1,0,0) and the corresponding quantities obtained from 
the oblique orientation. 

Our BCA code follows standard implementations^ It 
allows for arbitrary positions of the bulk atoms and is 
suitable for studying defect accumulation during multiple 
impact (although this feature is not used in the present 
work). Simple and well-tested forms of the screened 
Coulomb potential 26 and the inelastic processes22*2I have 
been chosen. All model parameters of the algorithm have 
been adjusted to Cu projectiles of a few keV, hitting a 
Cu single crystal^ These choices allowed for easy cali- 
bration and comparisons of our implementation against 
literature results. It should be emphasized, however, 
that the main focus of the present work is on generic 
results, which are of relevance for pattern formation of 
ion-sputtered surfaces of metals. 
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IV. 



SIMULATION RESULTS 



FIG. 1: Sample cascade originating from an impact of a 5 
keV Cu ion on a Cu crystal. The angle of incidence is 60°. 
The cube shown acts just as scale and has size 2.65 nm 3 , while 
the full lattice simulated has size (10.6 nm) 2 x 18 nm. 

We have performed BCA simulations of single ion im- 
pact on a plane Cu-surface at normal incidence with ve- 
locity v = - \v \e z . 

All the statistical information was obtained from en- 
sembles of 3000-6000 ion impacts per ensemble, which 
we generated for a single initial condition of the ion. The 
positions of the atoms making up the undisturbed solid 
were displaced from ideal lattice sites of a Cu single crys- 
tal (170000 atoms of an fee structure with lattice con- 
stant 3.61 A corresponding to a solid of 10.6 nm x 10.6 
nm and a bulk depth of 18 nm) by uncorrelated, Gaussian 
distributed displacements to account for thermal fluctu- 
ations. For each ion of an ensemble, an additional homo- 
geneous lateral random displacement was added, which 
was taken to be uniformly distributed within a square 
of edge length 1 lattice constant. Thus within every en- 
semble the ion hits upon macroscopically identical but 
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FIG. 2: Spatial distribution of ejected Cu atoms emerging 
from 6000 independent trials of hitting the (x, y) crystal sur- 
face (oriented in (1,0,0) direction) with a single 5 keV Cu 
ion at normal incidence. Distances are measured in units of 
a = 3.61 A. 

Fig- HI shows the surface distribution of all the ejected 
particles within an ensemble of 6000 cascades, each 
emerging from one incident 5 keV Cu ion (normal in- 
cidence) for the crystal in (1,0,0) orientation. Clearly, a 
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"hole" around the location of impact is visible. This is in 
contrast to what can be expected when applying Eq. . 
The scattering is almost rotational invariant, a slight 90° 
rational-invariant structure is visible, reflecting the lat- 
tice structure. To check whether the result is an artefact 
of the crystal orientation, we studied also a (58,32,39) 
surface, see Fig.|3J Although the plot exhibits slightly less 
structures, again only few particles are ejected near the 
point of penetration. Hence, we decided to concentrate 
on the oblique (58,32,39) orientation, because we want 
to study generic results irrespective of specific crystal ori- 
entations. In any case, for both surface orientations and 
small values of A, the angular average of data obtained 
from 



ri2d(p) 



1 



271-pA 



P+A 



pdp n 2 di.P,4>), (8) 



where p = \p\, turns out to be nearly undistinguishable. 
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FIG. 4: Probability of ejected particles vs distance p from 
point of ion incidence (measured in units of a = 3.61 A) de- 
termined from the data of Fig. . 
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FIG. 3: Spatial distribution of ejected Cu atoms emerging 
from 6000 independent trials of hitting the (a;, y) crystal sur- 
face (oriented in (58, 32, 39) direction) with a single 5 keV Cu 
ion at normal incidence. Distances are measured in units of 
a = 3.61 A. 
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FIG. 5: Surface density of mean energy of sputtered Cu 
atoms vs distance p (measured in units of a = 3.61 A) from 
point of ion incidence for 5 keV Cu ions on semi- log scale. 
The solid line is the best fit of the data to an exponential 
with a polynomial prefactor and corresponds to 0.297(p 2 — 
0.392p) exp( — 1.27p). The dotted line, which corresponds to 
a fit to a Gaussian, is obviously inadequate. 



Fig.0|shows the corresponding probability density p(p) 
per surface unit, averaged over all angles, of finding an 
ejected particle a distance p from the point of incidence of 
the Cu ion. This figure shows that the assumption on the 
ejection probability being distributed following a Gaus- 
sian distribution, hence leading to a maximum at p — 0, 
is not justified in case of crystals. This is in contrast 
to amorphous materials, where a more Gaussian- like dis- 
tribution of the ejection probability has been observed^, 
when using the simulation packet SRIM with the same 
ion/bulk parameters as above. 



Fig. [3] displays the corresponding angular average of 
the surface density e2d(p) of the energy of sputtered par- 
ticles. In this figure, we have also shown two Marquardt- 
Levenberg fits (f s — (ap 2 + bp) exp[— cp s ] with s = 1 
and s = 2) to the data. One can see that the decay of 
the energy density is not in accordance with a Gaussian, 
even when including a decay towards the point p = of 
penetration, as suggested by Eq. (JTJ. The data can be fit- 
ted well to a exponential decay with a simple polynomial 
prefactor. 

In the previous figure, we have studied the mean, hence 
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let us now consider characteristic features of the prob- 
ability density p(p, e) . In Fig. HJ1 we show the surface- 
integrated probability density p(e) = J d 2 p p(e,p). The 
behavior is remarkable close to a simple power law 



(b + e)<* 



(9) 



outside of a region of small e. Our best fit corresponds 
to a = 5.26, b = 5.03, a = 1.87. 
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FIG. 6: Probability density of energy, which is transported 
to the surface by a single collision cascade emerging from a 5 
keV Cu ion on a log-log scale. The solid line corresponds to 
5.259(5.035 + e) -1 ' 874 , which is the best fit to a simple power 
lawa/(6 + e) a . 

Fig. [7| displays the conditional probability density 
p{t\p) of energy at hxed p for different values of p. Sur- 
prisingly, the conditional density does not depend on p 
significantly. This shows that 



p(e,p) « p{p)p{e) oc 



an 2 d{p) 
(b + e) 2 



(10) 



outside of a region of very small distances. An immediate 
implication is that 



e2d(p) oc n 2 d(p) 



(11) 



so that the number of ejected particles and the energy 
deposited at the surface are proportional to each other, as 
assumed in the BH theory. However, another important 
implication is that the amount of energy transported to 
the surface is subjected to strong internal noise, which 
may limit the applicability of deterministic continuum 
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FIG. 7: Conditional probability density p(e\p) to find 
energy e of ejected particles keeping the distance p 
fixed. Different symbols correspond to different values 
p = 1.759, 2.914, 4.104, 5.277, 6.450, 7.623, 8.795, 9.968, 11.141 
As the data for different p are almost undistinguishable, dif- 
ferences not being significant, we need not provide a detailed 
legend. 
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FIG. 8: Distribution of ejected particles with different di- 
rections of velocity. Here k denotes the angle between the 
projection of particle velocity onto the surface and the vector 
between point of ion impact and point of particle ejection. 
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theories based on the average energy at the surface. This 
problem will be pursued elsewhere. 

To gain further insight into the structure of the en- 
semble of cascade particles at the surface, we have also 
studied correlations between the position p and the pro- 
jection of velocity, of sputtered particles onto the surface, 
Vsurf = v — (v ■ e z )e z . Fig.[S]shows the distribution of 
the angle n between p — po and v sur f, where po is the 
point at which the ion hits the surface and p is the po- 
sition of a cascade particle arriving at the surface with 
velocity v. The figures shows that, for typical collision 
cascades, most of the ejected particles move away from 
the point of ion impact. 



V. CONTINUUM APPROXIMATION TO 
ENERGY DEPOSITION 

Within Sigmund's approximation^ the rate at which 
the target is being eroded at an arbitrary point on the 
surface, is proportional to the total amount of energy de- 
posited there from ion collisions. In his theory for amor- 
phous or polycrystalline targets, an accurate description 
of the sputtering phenomena can be achieved by assum- 
ing that energy is deposited following the Gaussian dis- 
tribution QJ. 

Bradley and Harper (BH) 7 later employed this energy 
distribution in order to compute the local erosion veloc- 
ity at an arbitrary surface point O, allowing for gentle 
surface undulations. To perform the calculation, a new 
local reference frame is taken in which the z' axis is taken 
along the surface normal at O. The principal curvatures 
are assumed along the x 7 and y' axes, that are defined, re- 
spectively, as the direction orthogonal to z' that is in the 
plane defined by this axis and the ion trajectory and the 
remaining direction in order to make up a right-handed 
reference frame. Assuming that the radii of curvature 
at O, R x and R y , are much larger than the penetra- 
tion depth a, the surface height can be approximated to 
z'(x', y') = — + r - )' 111 or der to obtain the erosion 
velocity, we have to add up the total energy deposited at 
O from ions entering the whole target, expressing the ion 
flux and energy distribution in the latter reference frame, 
which is related with the one implicit in Q as 



x = x' cos(7o) + z' sin (70), 



y = y', 



(12) 



z'cos(7o) - x' sin (70), 



with 70 being the incidence angle formed between the ion 
trajectories and the surface normal at O. Accounting up 
to curvature corrections, the ion flux reads $(x',y') — 
<£>o ' (cos 70 — ■§- sin7o), where <&o is the constant nominal 
ion flux. Taking all this into account, the erosion velocity 
at O reads, finally, 

r+00 r+00 

v {jo,Rx,R y ) = A / / $(x',y')e s {x',y')dx' dy', 



where A is a proportionality constant relating deposited 
energy with the number of sputtered atoms, and the in- 
tegration limits are taken to infinity thanks to the fast 
decay of the energy distribution e s — eid- By expand- 
ing (| 1 31) to lowest non-trivial order in a/R x , a/R y <C 1, 
Bradley and Harper obtained 7 - 



v = A s Ae$ e" 



ro(7o) 



IM70) + r y(7o) 

Rx 



Ru 



(14) 



where To (70), r x (7o), and r y (7o) are functions that de- 
pend on the incidence angle 70, but also on features of 
the energy distribution such as a, a, and /3. 

Formula (|14|) enables computation of various relevant 
observables. Thus, the sputtering yield, ^(70), defined 
as the total number of sputtered atoms per incident 
ion, is easily related to vo by geometry as ^(70) = 
nvo{"fo)/(^o COS70), where n is the number of atoms per 
unit volume in the target. Assuming a planar interface, 
that is, in the R x , R y — ► 00, one is left with 



Y(7o) 



nvo(jo, R x — > 00, R y — > 00) 
$0 cos 70 



(15) 



Working with Sigmund's distribution, BH found- that 
y(7o) increases monotonously as a function of the in- 
digence angle 70, such that the maximum efficiency for 
erosion is achieved at grazing incidence, contrary to ex- 
perimental evidence for amorphous, polycrystalline, and 
crystalline targetsi^SLS This feature of BH's theory 
originates in a property of Sigmund's distribution l(T|l. 
whose maximum for deposition, r = (0, 0, —a), is located 
right at the surface under grazing incidence conditions. 
However, as is well known, there usually exists a value 
of 70 < 90° for which the yield is maximum, such that 
the sputtering efficiency decreases for larger angles of in- 
cidence, due to ions being reflected at the surface, an 
effect which is beyond Sigmund's approximations. 

Additional predictions on the morphology of the 
eroded target can be derived from l|14|) . Thus, r x (7o) 
is negative 7 - for small angles of incidence, which implies 
that the erosion velocity is larger at troughs (R x < 0) 
than at peaks (R x > 0), inducing a morphological insta- 
bility. Additional surface relaxation mechanisms exist, 
such as surface diffusion, that counteract this instability. 
Competition between the two opposing phenomena in- 
duce the emergence of a typical length scale, associated 
with the wavelength of the periodic ripple structure ap- 
pearing (for details, see Sec. IVIlj l. For 70 large enough, 
BH get T x > 0, while T y is always negative. Given that 
at small angles T x < T y , one obtains that the ripple 
crests are oriented perpendicular to the x' direction for 
incidences close to normal, whereas they are oriented per- 
pendicular to the y' direction for incidence angles larger 
than a critical one, 70 > 7 C such that T y < T x . Many 
experiments^^ have verified the validity of BH's theory 



(13) to describe ripple wavelength and orientation. 
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VI. MODIFIED ENERGY DISTRIBUTION 
FUNCTIONS 

The results of computer simulations within the BCA 
approximation, obtained in the previous sections for Cu 
ion bombardment of a Cu target are described by an en- 
ergy distribution that differs substantially from that ob- 
tained by Sigmund in the case of polycrystalline or amor- 
phous substrates. Using cylindrical coordinates around 
the ion trajectory, as in previous sections, we have [recall 
Fig. and Eq. JT^ above] 



e e (p, z) = N e e(p 2 + cp)e a *v e 



(16) 



where N e = [(2n) 3 / 2 a z (6a xy + 2ca xy )]~ 1 is a normaliza- 
tion constant. Values for c and o~ xy that best fit sim- 
ulation results were c = —0.392, a xy — 0.787, see Fig. 
El Note two main differences of distribution (fT?)|l to Sig- 
mund's distribution QJ: decay here is slower [exponential 
as compared to Gaussian, thus the subscripts in HIGH ] in 
the plane perpendicular to the ion trajectory, and en- 
ergy deposition is null along the ion trajectory itself. On 
the other hand, distribution H16|) is unphysical. since 
c < leads to negative probabilities for small p values. 
For this reason, and in order to facilitate analytical re- 
sults [unavailable for 1(16(1 in the physical case of a two 
dimensional target, see below], we will also consider the 
following modified Gaussian (thus the subscripts) distri- 
bution 



e g (p,z) = N g ep 2 e 2 °*y e 
N g = [2(2^) 3 / 2 a xy a z \ 



(17) 



that shares with 1(16(1 inducing zero energy deposition 
along the ion trajectory p — 0, but is otherwise Gaus- 
sian in all three directions far enough from the ion path. 



A. One-dimensional interfaces 

In order to develop intuition about morphological pre- 
dictions from 1(16(1 . 1(17(1 . that are based on analytical re- 
sults, we consider first the (non-physical) case of a one- 
dimensional target, whose surface height is described by a 
single variable function z'(x'). These results will be then 
compared to the analogous ones by Bradley and Harper, 
which will allow us to assess differences due to the new 
form of the energy distribution, mostly to the fact that 
in our case no energy is deposited along ion trajectories. 

For a one-dimensional target, distribution 1(17(1 reads 



e ld (x,z) = N g ex 2 e 



N ld = (2na z a 



3\-l 



(18) 



Writing the local velocity of erosion in terms of Tp' 1 ^ 
rg' ld analogous to Eq. lfH)l . we obtain 



vo 



N} d Ae<S> e~ 



(19) 



where the full expressions for Tp' 1 " 1 and T 9 x ,ld as functions 
of 7q, a, a xi and a z can be found in Appendix I A II 

On the other hand, distribution 1(16(1 reads, for a one- 
dimensional interface, 



e\ d (x,z) = N e e{x 2 + c\x\)e °* e 

2M-1 



(20) 



N. 



Id 



In this case, the prediction for the local velocity of erosion 
has a shape that is similar to (|18fl . albeit with more com- 
plex coefficients, whose detailed analytical expressions 
are again left to Appendix I A 21 



v = AT ld Ae$ e 



-,e,ld 



Ad 



R x 



(21) 



In Fig. El we plot the normalized (to the corresponding 
values for normal incidence) sputtering yields ^(70) ob- 
tained through (|15f) for the modified Gaussian (|19fl and 
exponential (|21|l distributions. For the sake of reference, 
the BH yield is also shown. We can see that for both 
modified distributions, i|18|) and H20|) . the corresponding 
yields feature maxima before grazing incidence, as a dif- 
ference to the BH curve. This is in agreement with exper- 
imental data^i^S and is due to the fact that maxima 
of energy deposition are not along the ion trajectory for 
these distributions but, rather, at a certain finite dis- 
tance from it, what makes grazing incidence not be the 
most efficient for sputtering. The fact (seen in Fig. EJ) 
that the yield is negative for large incidence angles 70, 
as computed using the exponential distribution, is due 
to Eq. (|2*U|) taking negative values for small distances to 
the ion path. As will be seen below, this is an artifact 
of the one-dimensional approximation, as is the fact that 
the yield computed from the modified Gaussian distribu- 
tion l|18|) vanishes for 70 = 90°. From the figure we can 
also conclude that qualitative behaviors of distributions 
(|18fl and (|20fl are similar, the advantage of the first one 
being its greater analytical simplicity. Parameters em- 
ployed in Fig. El are typical for Cu ion bombardment of 
Cu for energies in the range of a few keV, as confirmed 
by TRIM/SPJM simulations^ 

B. Two-dimensional interfaces 

Naturally, the physically relevant case is bombardment 
of two-dimensional targets. In this case, the analysis is 
more complex, to the extreme that no closed expressions 
analogous of those found previously for the exponential 
distribution 1)16(1 that best fits our BCA simulation data. 
Results for this distribution will be provided from nu- 
merical solutions of (|13|) using (|16|) . On the other hand, 
we have seen in the Id case that distributions 1(18(1 and 
1(20(1 lead to similar qualitative results. For the 2d case, 
expression (fH3|> using (fTTjl leads to closed analytical ex- 
pressions for the coefficients in 
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FIG. 9: Normalized sputtering yield y(7o)/V(0) as a func- 
tion of incidence angle 70, for the various one-dimensional en- 
ergy distributions. Dashed line: Bradley-Harper for a — 3.8 
nm, a = 2.2 nm, (3 — 1.5 nm. Solid line: modified Gaussian, 
Eq. CHJi f° r a = 3.8 nm, a z — 2.2 nm, a x = 1.5 nm. Dotted 
line: exponential, Eq. (I20L for a — 3.8 nm, a z — 2.2 nm, 
era, = 0.787 nm, c = -0.392 nm. 
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FIG. 10: Normalized sputtering yield Y(-yo)/Y(0) as a func- 
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energy distributions. Dashed line: Bradley-Harper, Eq. 
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nm. Dotted line: exponential, Eq. 11611 . for a = 3.8 nm, 
a z = 2.2 nm, a x = 0.787 nm, c = -0.392 nm. 



that can be found in Appendix IA 31 In Fig. EH we see 

again that the sputtering yield for both modified distri- 
butions (|16fl and l|17|l have maxima for incidence angles 
clearly smaller than grazing. Moreover, the yields are 
positive and non-zero for all values of 70, and amount to 
large sputtering rates, as found in experimentsi2ii22ii£. 
In the present two-dimensional case, for grazing incidence 
the radial component of the energy distribution vanishes 
at the point of impact with the surface, but not at finite 
distances from it, which implies that after surface inte- 
gration the total deposited energy is non-zero and the 
yield is positive. Again, Fig.lTfjIshows similar qualitative 
behaviors for both modified distributions, the modified 
Gaussian having the advantage of leading to closed ana- 
lytical results. For normal incidence, the x <-> y symme- 
try is restored, and thus ^(70) and T y (70) must coincide. 
From Eq. we obtain r»(0) = rg(Q) = -4jraa% y /a*, 
whose coincidence with the numerical results for distri- 
bution i|16|) has been confirmed. In Figs.^J^l and 1131 
we present results for the "surface tension" coefficients 
T x and T y for the various two-dimensional distributions, 
that are normalized by their corresponding absolute val- 
ues for normal incidence 70 = 0. We see in Figs. fTTH IT^l 
that T x is smaller than T y for incidence angles 70 < 7 C 
and that T y is always negative, similarly to the BH case 
(Fig. One of the successes of BH's theory lies in 

its description of the orientation of the ripple structure 
for different ion incidence angles 70. Here we see that, al- 
though distributions 116|) and (|17|l . lead to quite different 
sputtering yields as compared to Sigmund's distribution, 
the qualitative behavior of coefficients T x and T y is quite 
similar to that found by BH. Since experimental results 
are in good agreement with BH for metals, the modified 
Gaussian distribution (|17|) seems a good choice for the 



corresponding analytical description, with the advantage 
over the exponential distribution i|16|) of being physically 
sound for small distances to the point of penetrationii^ 

30 ! 
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FIG. 11: Normalized values of T x and T y for the distribution 
Q using the same parameter values as in Fig. 1101 



VII. CONTINUUM EQUATION FOR THE 
SURFACE HEIGHT 

Following the pioneering approach by Bradley and 
Harper, we can derive an evolution (differential) equation 
for the surface height, starting from the equation for the 
erosion velocity. We consider a laboratory frame of refer- 
ence (X, Y, Z), defined as follows: the Z axis is chosen 
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FIG. 12: Normalized values of T% and F y for the distribution 
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FIG. 13: Normalized values of T% and Ty for the distribution 
1161 1 using the same parameter values as in Fig. 1101 



to be normal to the initial planar surface. The incoming 
beam direction forms an angle 9 with Z, and both direc- 
tion define a plane where the X axis lies. Finally, the 
Y axis is perpendicular to the X and Z directions. We 
describe by h(X,Y,t) the surface height at time t above 
point (A, Y) on reference plane of the unbombarded sub- 
strate, and assume that it varies slowly enough so we can 
work to first order in the derivates. In this way we may 
approximated 7 o = 6-%, £ = -0, £ = -0. 
The velocity of erosion of the surface height h is provided 
by the erosion rate vo, and we thus get: 



ar o (0) dh 



T 



d 2 h 



d 2 h 



89 dX x dX 2 v dY 2 ' 



(23) 



where in our normalization F is a proportionality con- 
stant between vo and r , T x , T y , that can be found in the 
Appendix. Considering a periodic perturbation to the 
planar surface h(X, Y,t=Q) = Ae^ klX+k2Y \ and substi- 



tuting this expression into Eq. I|23[) . the surface profile 
evolves as 



h(X,Y,t) = -T t + Ae rt e l{klX+k2Y - ut \ 
r = — r^fcj — Tyk 2 , 
lo = —T ki. 



(24) 



If T x and/or T y are negative, there will be values for 
the wave- vector (ki,k,2) of the perturbation that make 
it grow exponentially. This behavior is a reflection 
of the well-known physical instability leading to ripple 
format ion j£>2P. due to the curvature dependence of the 
erosion velocity, that is larger in surface troughs than 
in surface protrusions. The observed ripple wavelength 
arises when additional smoothing mcchanims such as sur- 
face diffusion exist that compete with the sputter insta- 
bility, leading to selection of a specific length-scale. Tak- 
ing these mechanisms into account Eq. I|23|) reads 



dh 



— = F 



-r (6) 



dT (6) dh 



d 2 h 



89 dX x dX 2 
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d 2 h 



V dY 2 

(25) 



where, in principle, B is a thermally activated coeffi- 
cient which depends on the surface self-diffusivity D s , 
the free energy per unit area 7 and the number of 
atoms per unit area moving across the surface a as 
B = 2D s ja/{n 2 k B T). In this case, r = -T x k 2 - Y y k\ - 
B{k 2 + k 2 ) 2 , and there is only a band of unstable per- 
turbations. The observed ripple wavelength £ is provided 
by the wave-vector which has the largest positive value 

of r, and is proportional to or \J F\r wn en 

T x <T V < (9 < 9 C ) or T y <T X (9 > 9 C ), respectively. 

As we have seen in the previous section, the behavior of 
F x and T y is similar for all the cases considered, while the 
qualitative behavior of the yield is quite different to that 
found by Bradley and Harper for amorphous or policrys- 
tallinc substrates. Moreover, since dependences of the 

ripple wavelength I cx ^/-pj^y on parameters such as ion 
flux, $0, temperature, or average ion energy, e, are due to 
those in the constants F and B, and these are the same 
as those in BH (see Appendix), Eq. (|2*5|l predicts these 
for Cu to be (qualitatively) the same as obtained from 
BH theory^ Note this is also the case in the presence of 
non-thermal surface diffusion, in which, similarly to BH, 9 
the constant B has no dependence on temperature and 
is, rather, proportional to F. 



VIII. SUMMARY AND OUTLOOK 

We have studied numerically the sputtering process of 
Cu ions on Cu fee crystals by means of the binary col- 
lision approximation. We have analyzed the distribu- 
tion of sputtered particles and their energies, and found 
significant deviations from Sigmund's formula, which is 
traditionally employed to study the sputtering process 
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in the framework of continuum theories, as applied to 
amorphous and policrystalline substrates. In particular, 
we find that near the point where the ion penetrates 
the target, the sputter probability goes to zero, while 
the Bradley-Harper/Sigmund theory predicts maximum 
sputtering at that point. 

We have fitted heuristic functions to our data. We 
find that an exponential (rather than Gaussian as in Sig- 
mund's theory) decay with a combination of a quadratic 
and a linear prefactor fits the data well. The main phys- 
ical effect, namely, the "hole" near the point of penetra- 
tion, can be reproduced also qualitatively using a Gaus- 
sia distribution with with a quadratic prefactor, that is 
physically better defined than the heuristic fitting distri- 
bution, and lends itself to exact results. We have per- 
formed analytical calculations of the local erosion veloc- 
ity following the Bradley-Harper approach for one- and 
two-dimensional surfaces, for both types of modified dis- 
tributions (for the two-dimensional exponential distribu- 
tion, the equation could be solved only numerically) . We 
find that the sputter yield is qualitatively different as 
compared the the BH approach. As a function of the 
angle of incidence, the yield exhibits a maximum at an 
intermediate angle, and then decreases when approaching 
grazing incidence. This is in good agreement with exper- 
imental findings, in marked contrast with the analogous 
BH result using Sigmund's distribution, even without im- 
plementing explicitly reflection of the ions for grazing in- 
cidence, which is usually regarded as the main cause for 
the decay of the yield at grazing incidence. Finally, we 
have computed also the ripple orientation-determining 
parameters T x , T y , usually referred to in this context as 
effective surface tension parameters. These turn out to 
be only slightly modified with respect to the BH the- 
ory, and lead to a qualitatively similar pattern formation 
process. Dependencies of the ripple wavelength on phe- 
nomenological parameters, such as ion flux, ion average 
energy, and temperature are as in BH theory. 9 Since the 
influence of non-linearities on ripple characteristics is still 
under debate even within Sigmund's theory proper, we 
have not considered this type of effects here. At any 
rate, the same type of non-linear terms would appear in 
the interface Eq. I|25[) as compared to the corresponding 
equation for amorphous or policrystalline substrates^S 

Thus, as a general conclusion on pattern formation 
by ion-beam sputtering, our results justify the similar- 
ities found in experiments on metals, to the analogous 
processes in amorphous or amorphizable materials, and 
point to potential quantitative differences that would 
possibly merit further studies. Additional features of rip- 
ple formation in metals, such as their existence for normal 
incidence or change of orientation with temperature 3,15,16 
are not explained by the special properties of the colli- 
sion cascades in these systems that we have studied here 
but, rather, by the special properties of surface diffusion 
in such anisotropic substrates. 

Regarding future work, it would be also interesting to 
see whether the hole near the point of penetration can 



be found in experiments, and/or in more detailed sim- 
ulations (such as e.g. by Molecular Dynamics). To our 
knowledge, no analysis of single-ion impacts on metals 
exist so far. Furthermore, it would be worth incorporat- 
ing the modified energy distribution into existing simple 
Monte Carlo models of surface sputtering, such as those 
m 

Refs.ElE 

in order to improve their description of 
erosion processes in metallic substrates, specially at the 
large distance and long time regime for which this type 
of models is particularly suited. 
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APPENDIX A: ANALYTIC EXPRESSIONS FOR 
COEFFICIENTS IN THE EROSION VELOCITY 

In this appendix, we provide the full expressions for 
the coefficients appearing in various expressions for the 
surface erosion velocity, Eqs. I|21|) . (|19l) . and l|22|) . that 
have been computed analytically for those energy distri- 
butions for which such type of results are achievable. 
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2. One-dimensional exponential distribution 
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